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Abstract. Agents' heterogeneity is recognized as a driver mechanism for the 
persistence of financial volatility. We focus on the multiplicity of investment 
strategies' horizons, we embed this concept in a continuous time stochastic 
volatility framework and prove that a parsimonious, two-scale version effectively 
captures the long memory as measured from the real data. Since estimating 
parameters in a stochastic volatility model is challenging, we introduce a robust 
methodology based on the Generalized Method of Moments supported by a 
heuristic selection of the orthogonal conditions. In addition to the volatility 
clustering, the estimated model also captures other relevant stylized facts, 
emerging as a minimal but realistic and complete framework for modelling 
financial time series. 
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1. Introduction 

In 1963 [TJ [5] Benoit Mandelbrot refers to the volatility clustering as "large changes 
tend to be followed by large changes, of either sign, and small changes tend to be 
followed by small changes". Since then this effect has remained one of the most 
intriguing properties exhibited by financial time series. In the early Nineties the long 
memory property of absolute stock market returns was independently investigated 
by [3] and [J]. In the former work, after amending absolute price changes from 
the heteroscedasticity due to seasonal effects, the authors find a persistent positive 
autocorrelation declining hyperbolically with the time lag. In the latter, analysing 
the daily closing prices of Standard&Poor 500 index for the time span January 3 1928 
- August 30 1991, Ding and collaborators study the power correlation of absolute 
returns \r t \ d for positive d, finding a strong persistence especially for d close to one. 

The slow decay of the volatility can be ascribed to two rather different 
mechanisms. Agent Based Models provide a first explanatory framework, where 
macroscopic evidences are explained in terms of microscopic interactions among 
market participants. As clarified in the seminal papers |6] the alternation of 
the economic agents between chartist and fundamentalist regime can be identified as 
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the source of the observed volatility clustering, an empirical signature of persistence. 
The same mechanism leading to the previous regime switching is further investigated 
in pj [5] , where the minimal assumptions required for an agent based model to capture 
the empirical stylized facts are identified. In a different approach [9 persistence is 
induced by the coexistence of agents differing in their perceptions of the market, 
risk profiles, institutional constraints, degree of information, prior beliefs, and other 
characteristics such as geographical locations. In [10] the role of heterogeneous time 
horizons for the investment strategies is specifically addressed. In [TT] the daily, weekly 
and monthly time scales are isolated as the relevant ones, while the first direct evidence 
of these three scales as well as an attempt to capture them with an ARCH model is 
provided in [12]. As a major achievement of the latter work we see that a small subset 
of time scales succeeds in capturing the long run behaviour of the squared return 
correlation. Interestingly, those horizons reflect typical time scales of the human 
activity, which noticeably follow a pseudo-geometric progression [13.. Generalizing 
the concept of a finite mixture of time scales to a continuum of agents, an attractive 
intuition is that the integrated effect of exponential heterogeneous strategies may lead 
to persistence. On a formal basis, this amounts to expressing the correlation function 
as 

/•l/r mJn 

C(t) = / exp(-r/r agcnt )p(l/r agent )d(l/r agcnt ) , 

Jo 

which at the leading order for r — > +oo is determined by the behaviour of the density 
p(l/r ag ent) around the origin. Indeed, by virtue of Watson's Lemma, we obtain 
C(t) ~ l/r 1+a provided that p(l/r agont ) ~ T~* nt with a > -1. 

As far as the distributional properties of the volatility proxies are concerned, 
in [14] the inverse gamma distribution is identified as an effective approximation 
for both the low and high volatility regimes. The simplest model reproducing this 
distribution as a result of a volatility feedback effect corresponds to an ARCH-likc 
equation which, in the continuous time limit, reads as a Langevin equation 

da 

= -4°" - a oo) + v ° Cw , 

with k, Coo, rj positive constants. For this specific case, the stationary distribution of 
the volatility has the form of an inverse gamma 

Y e- A / ff 
I» <J X + V 

with v = 1 + 2k/ i] 2 and A = 2Ka 00 /rj 2 . In the following Section we propose an 
approach inspired by this evidence about the volatility, as well as by the idea of a 
mixture of heterogeneous investment horizons, in a spirit similar to the Heston multi- 
factor model [15] , 

The remainder of the paper is organized as follows: in Section [3] we derive 
analytical expressions for the leverage and volatility autocorrelation, while in Section[4] 
we detail a calibration procedure which is inspired by the Generalized Method of 
Moments. We conclude in Section [5] The analytical derivations are postponed in the 
three Appendices at the end of the paper. 
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2. The model 

A quite general expression for the asset price at time t, reminiscent of the Geometric 
Brownian motion paradigm, is given by 

S t = So exp (fit + X t ) 

where X t is the stochastic centred log-return and /i a constant drift coefficient. In [16] 
we assume that the time evolution of X t can be modelled in terms of the stochastic 
differential equation (SDE) 

dX t = a t dW t x , (1) 

where o~ t is the instantaneous volatility of the price and dVF t x the increment of a 
standard Wiener process. Since Xq is equal to zero, we also have E [X t ] = and 
E [In St — In So] = [it for any t. A common choice accounting for the stochastic 
behaviour of the volatility, as measured by suitable proxies, is cr f = o~(Y t ) as a function 
of an unobserved driving process Y t . General financial considerations regarding the 
mean-reverting behaviour of the volatility process lead to a second SDE of the form 

dY t = ~K Y {Y t - Voo ) At + y/X(Y t ) dWf , (2) 

with Ky = 1/ry- > 0, and y^ > 0. In [TB] T,(Y t ) is equal to o Y Y^ with oy > 0, 
from which it follows that a t is proportional to Y t . This choice leads to an inverse 
gamma stationary distribution with shape and scale parameters v = 1 + 2ny /o~ Y 
and A = 2 Ky y 'oo / 'o~y, respectively; in light of the considerations presented in the 
Introduction, this is dictated by the will of recovering the most effective statistical 
description of the volatility distribution. Different choices for £ have been suggested in 
the literature and among the most popular ones it is worth mentioning the Heston |17j 
and Stein-Stein [TB] models. For a complete overview of continuous time models 
as well as widely employed discrete time approaches like ARCH, GARCH and their 
generalizations, we suggest the handbook about financial time series |19j . 

Following the spirit of the Introduction, in this paper we extend the model given 
by (1) and (2) allowing the instantaneous volatility a t to depend on multiple stochastic 
unobserved factors. Before stating explicitly our model's equations, it is worth noting 
that such a generalization is inspired by the multi-factor stochastic volatility model 
introduced by Bates in [TS] as a possible description of the S&P 500 futures price, 
and later revisited in [20] as a model for the dynamics of the volatility smirk in the 
option pricing context. In the generalized form used in [21] , the multi-factor model 
with jumps reads 

N 

dX t =^2^dWl + dJ t x 

i=l 

dYi = - Ki (Y? - yl) dt + m ^dWl +N +dJj, i = 1, . . . , N 

(3) 

where Wt,...,Wf N is a multivariate possibly correlated Brownian motion and 
{J x , J 1 , . . . , J N } is a multivariate possibly correlated Poisson process with constant 
intensities. In principle, each factor may be linked to the sensitivity of the economic 
agents to different investment horizons, and in light of this heterogeneity the modelling 
could reflect N volatility components. Our starting point is a special case of the 
dynamics (3): we take N — 2, discard the jump contribution and assume perfect 
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correlation between the Wiener processes of the two factors, corr(Wf,Wf) = 1. 
However, at variance with equation (3) and as a major contribution of our paper, 
we consider inverse gamma driving factors, each one being described by the same 
mean-reverting dynamics provided in (2) with S proportional to the squared process. 
Ultimately the model we are going to analyse reduces to 

dX t = Y t dW t x + Z t dW t x 



dY t = -K Y (Y t - Voo) dt + cry Y t dW t 



dZ t = -K Z (Z t - zoo) dt + g z Z t dWf , (4) 

where we impose the initial time conditions X t= o = Xq = 0, Y t= t — yo > and 
Zt=t = z o > 0, with k y — 1/ty > 0, and k z = 1/tz > 0. We also indicate 
Vy = 1 + 2k y /o~y an d v z = 1 + %kz/o~z the tail exponents of the inverse gamma 
stationary distributions of Y t and Z t . The correlation structure among the three 
Brownian motions is described by the following matrix 




It has to be noted that we assume different starting times for the volatility factors 
and the return process, according to what is done in [16 . Indeed, as the processes Y t 
and Z t are unobserved factors and we are mainly concerned with their dynamics at 
the stationary state, we assume they start at t < in the past and we recover the 
stationary limit by letting to —> — oo. On the other hand, X t represents the observed 
(detrended) logarithmic increment of the price for a fixed time lag and, therefore, it 
seems natural to take the spot time t = as a starting time for this lag. 

Some considerations are due regarding our choice of the factors specification. 
In [TB] the single factor Y t corresponds (up to a constant) to the instantaneous volatility 
itself. As such, Y t has a clear interpretation and its dynamics is chosen specifically 
with the intent to accommodate the distributional properties of the volatility observed 
in the reality. Here, in the spirit of the factor model ([3]), the evolution of log-returns is 
given in terms of two additive factors; it follows that Y t and Z t can not be interpreted, 
separately, as the return volatility (or possibly the variance as in ARCH/GARCH 
models), but as the underlying unobserved factors. The choice of full correlation 
between W} and Wf allows us to introduce formally a t = Y t + Z t and to motivate the 
inverse gamma dynamics which drives the factors. Implicitly, this means that we give 
up recovering exactly an asymptotic inverse gamma law for at and we give priority 
to capturing the observed, long range memory of the squared return correlation. 
Nevertheless, we expect the tail asymptotic to be preserved under suitable assumptions 
(see discussion at the end of this Section). Finally, we observe that generalization with 
more that two factors is straightforward, but cumbersome, and would greatly simplify 
if we assume corr(W/, W?) — in Q. However, the specific purpose of this work is 
to show how our minimal choice is indeed able to capture the very consequences of 
heterogeneity. 

From |16j we know that a negative pxy suffices to accommodate the observed 
short range scaling of the return-volatility correlation; in |Appcndix B| and in the 
numerical section we set pxz equal to zero to prevent Z t from impacting the leverage. 
Nonetheless, in what follows we derive the relation between the factors behaviour 
and the moments of X t under the general case of non trivial correlations between 
the Brownian motions. As we show in Appendix A the structure of the model Q 
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allows to compute the moments of the probability density function (PDF) of X t at 
all times t recursively. After cumbersome calculations, and by exploiting Ito's Lemma 
to compute the cross correlations between the two volatility factors, it can be verified 
that the moments of X can be expressed always as a superposition of exponential 
functions of (f — to) 



where the constants read F m ^ n 
m(m-l)o-y/2, and F% = 



E [A™] = Bijft J*' z o) ex P ( F iAt ~ to)) , 

t,j=0; i-\-j<n 

F m + F n + mnpYz tJv'y" z ■ 
n(n 



(5) 



2 crL, with Fl 



-Kym 



l)o~ 2 z /2. The coefficients depend on the 



time lag f; more precisely, due to the linearity of the ODEs (A.2), they correspond to 



a combination of exponential terms weighted by polynomials in t. 

In the following, we report the explicit expressions of the coefficients H\ n ) (t; yo, z ) 



for the case n = 2 (the constants k^' n \ which depend on the initial conditions yo, z , 



are defined recursively in Appendix A|) 
r(2) 



H o,a - 

^1,0 - 



rr(2) _ 
fn i — 



1.(2,0) 
"■1,0 



v 0,0 



2^1,6 "* 



,(2,0) 



•2fe 



(i,i) 



1 - exp (-Fi fi t) 
Fi,o 

1 - exp (-Fo,it) 



(2) _ 



jj(2) 
"1,1 



(2,0) 
2,0 



1 - exp (-F 2 fit) 



Fo 



H. 



(2) _ 
0,2 " 



, (i,i) l- exp (- 

exp (- 
-Fb,2 

r(") 



1.(0,2) 



1 



^0,2*) 



Since t is finite, the coefficients H\J are finite quantities themselves, and all the 
relevant information about the behaviour of E [A"] in the stationary limit of Y and Z 
is retained by the to- ex P orLeil tials in Equation |5]). Given that F — 0, if all the F^ 
for i,j = 0, ... ,n with i + j < n are negative, E [X™] is finite in the stationary limit 
to — > — oo, otherwise it diverges indicating the emergence of fat tails in the PDF pt(x) 
of Xt ■ In the latter case the tail behaviour would be compatible with an hyperbolic 
scaling with a tail exponent smaller than the order of the lowest diverging moment. 

In |16j the hyperbolic scaling of pt(x) is induced by the power-law tail of the 
asymptotic (inverse gamma) distribution of the volatility, and a simple relation exists 
between the tail exponent of the latter and the order of the first diverging moment 
of Pt(x). In the present case, the asymptotic distribution of at is that of the sum of 
the two factors Y and Z , both inverse gamma distributed with tail indices Vy and v z 
respectively. In the limit of Y independent of ^f] the distribution of the sum behaves 
as a power-law with tail index v a = min-jVy,^} = see e.g. [221 123] . Therefore, 
the same mechanism discussed in |16j . which triggers the divergence of the return 
moments, applies here asymptotically in the absolute value |A|, and the PDF pt(x) 
manifests a decay compatible with a power-law scaling with tail index determined by 
the value of v a . 

% In fact, we make this assumption when estimating the model from the empirical data in Section |4j 
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The theoretical PDF is also compliant with the more basic properties of the 
returns. In particular, the stationary limit of (pL and the expression of Hqq, show that 
the variance is linear in the return time lag t. furthermore, the explicit expressions of 
the functions flg o an< ^ d > n °t reported here for the sake of parsimony, would also 
reveal that the skewness and kurtosis vanish in the limit of large t, according to the 
observed Gaussian-like shape of the distribution for large time horizons. 

3. Non linear dependence 

In this Section we discuss the main properties of the return correlation structure 
predicted by model Q, focusing on the return- volatility and the squared-return 
correlation functions. 

Model Q inherits from the class of stochastic volatility models the important 
property of absence of serial correlation, which is verified empirically with good 
approximation. Despite this, financial returns can not just follow a random 
walk process, since this would imply independent and identically distributed price 
increments. Any non linear function of the returns would then exhibit zero 
autocorrelation, a property that simply does not hold in practice. Empirical evidences 
of this violation are the leverage effect and the volatility clustering. The former 
refers to the negative correlation between past returns and the future instantaneous 
volatility, measuring the tendency of the market volatility to increase after a price 
downfall [24, 22, 25j; volatility clustering is usually expressed in terms of the persistent 
correlation between squared returns or logarithm of absolute returns, implying that 
large variations are more likely to be followed by large than small ones [3, 26, 27, 28, 29 . 
For a survey of contributions on the same topics from the econometric community, we 
refer the interested readers to the reference list in jT3j [30] . Our model deals explicitly 
with these non linear correlation functions; their expressions, whose derivation we 
postpone to Appendix A| - |A"ppcndix~C| provide a valuable analytical characterization 



of model |4). This information can be exploited for the calibration of the model 
from empirical data, but it is also crucial to grasp the relevant information about the 
process time scaling, as we discuss in the rest of this Section. 

3.1. Leverage effect 

The leverage, a measure of the correlation between returns and volatility, is usually 
defined as C(r;t) = E [dX t dX? +T ] /E [dXf] 2 . Empirically and for arbitrary i, 
£(r; t) has been found to be negative and exponentially decaying for positive r 
and approximately zero otherwise: a correlation exists between past returns and the 
volatility in the future and not vice versa. 

For our case, a finite time, exact expression is derived in |Appcndix B[ Equation 



(B.2| reveals that the leverage function is characterized by the superposition of three 
exponential functions, with different characteristic times T£, t^c, t<, which are 
ordered according to the following hierarchy: 

2 vy - 1 

2 — TyCTy Vy — 2 
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H- 


4) 


<41 






2 



2T T C = 

TCO-y 



Vy - 

2v Y 



If vy — > 3 + , converges to 2ry, while for j/y — > +oo we have that tc goes to Ty. The 
time scale t^£ is strictly smaller than tc', however, we are implicitly assuming that 
the characteristic time of the Z process accounts for the volatility persistence, that is 
tz>T£, implying that t^c is expected to be only slightly smaller than the leverage 
scale. Ultimately, if vy — > 3 + , t< converges to 2tc/3, while under the Gaussian limit 
we have that r< converges to tc/2. In fact, the three leverage scales are constrained in 
a narrow range, which empirically has been found to be of order ten days for indexes, 
or even larger for single stocks [22]. 



3.2. Autocorrelation function of squared increments 

The volatility clustering is commonly measured by the quantity E [dX t 2 dA 2 +T ] and 
the volatility autocorrelation can be estimated in terms of the following normalized 
quantity 

E [dXf dX? +T ] - E [dX?] E [dXf 1 
/Var[dA 2 ] Var[dX t 2 H 



A(t; t) = - ^ 7^ t+rJ = ^ J = ¥ U±lL . (6) 



To complete our analytical characterization of the two-factor model Q , in Appendix 



[C] we derive the explicit expression of ^. At variance with the one previously 
found in 1 6 an d which is unable to capture the persistence of the volatility, the 
expression (C.7| features five different exponential scales r^ -1 '" ' 5 '*. Similarly to the 
leverage, the characteristic times are organized in a hierarchy as follows 

(1) 1 TjC 

T A = ~^Z = ~ TZ = T > Z > TZ ' 
r 2 T Y 



1 



.(2) _ 
A — 

-(3) _ 



1 

w 



= Ty <T C , 



(4) 1 T Z 

Ta = -(f? + F z) = ^T^ Ty<Ty 

T (5) = _ 1 = IC 

A Fj 2 ' 



For vy varying in (4, +oo), Ty is inferiorly bounded by 2x^/3, while the upper bound 
is given by tc- Therefore r>^ ranges between tz and 3r^/2, and we can conclude that 
the previous five scales indeed cluster into two groups, a long-range and a short-range 
one: the first set is {Tz,Tyz}, whose typical scale is given by r^, while the second 
one contains the three remaining scales, superiorly bounded by T£ and of order ry. 
Ultimately, we can appreciate the very reason why model Q has been enriched by 
a second factor process Z t w.r.t. the one proposed and discussed in [TS]. Through 
the coupling provided by pxy ^ 0, Y t is entirely responsible for the emergence of the 
leverage; conversely, the Brownian motion driving Z t is decoupled from , can not 
interfere with the leverage, and can not constrain its hierarchy of time scales, as it 
happens in [16] . In model Q Z t provides the degree of freedom required to capture the 
persistence of volatility. It is not difficult to imagine that extra volatility factors would 
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induce a new plethora of time scales. However, even though the analytical tractability 
would be preserved, this would come at the cost of an overwhelming burden of messy 
calculations. 

4. Calibration via Generalized Method of Moments 

In this Section, we propose and discuss an application of the Generalized Method of 
Moments to the estimation of the model's parameters from en empirical time series of 
price increments. 

The stochastic model Q is characterized by twelve free parameters, ry, y^, yo, 
o-y, tz. Zoo, zq, o~ 2 z , pxy, Pxz, Pyz, and /i. Estimating parameters in a stochastic 
volatility model is a challenging task. This is primarily due to the latency of the 
volatility state variable. Indeed, in different approaches to volatility modelling, like 
ARCH and GARCH models, the likelihood function is readily available. This problem 
has inspired many scholars, and there is a specialized literature on computationally 
intensive methods mimicking likelihood-based inference. In general, these belong to 
the class of non linear filtering methods, and among possible approaches we mention 
Kalman filters, Particle filters, and Monte Carlo Markov Chain approaches. For 
more techniques and further discussion we refer the reader to the handbook [19,. 
Here we take the opportunity to quote the interesting proposal discussed in [3T] and 
rooted on the spectral approach to nonlinear filtering. A relatively simpler approach 
to estimation, which does not rely on any ad hoc approximation of the density of 
returns, is based on the computable moments of the model. For continuous-time 
stochastic volatility models, it is generally very hard to derive closed form solutions 
for the return moments, but this is not the case for the model under consideration. 
For this reason, we follow a methodology inspired by the Generalized Method of 
Moments (GMM). An introduction to the GMM, based on Hansen's formulation of 
the estimation problem [32], is provided by [33] in Chapter 14. Given T observations 
{W t } for t — 1,...,T, each one being an h dimensional vector, and a vector 
6 R k of unknown parameters, in order to apply GMM there should be a function 
h(0,W t ) : R k x R h -)• W characterized by the property that 

E[h(0,W t )] = O. (7) 

These r equalities are usually described as orthogonality conditions. The basic idea of 
GMM is to replace these conditions with sample averages and to solve the following 
optimization problem 

8 = aga|^h(8,W t )j fi" 1 (~ £ h(0, , (8) 

where £1? is a positive-definite weighting matrix depending on the available data set 
and on the value of 9 itself. The practical procedure is the one which follows. An 

initial estimate 6 is obtained by minimizing the previous quantity with an arbitrary 
choice of CIt, e -g- = Irxr- Supposing that h(0,W t ) is serially uncorrelated, the 

estimate (O) is then used in 

f2 T = if]h(0 (O) ,W t )h t (0 (O) ,W t ) 
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to arrive to a new GMM estimate 9 . This process can be iterated until an arbitrary 
stopping criterion is invoked [§J If 9 denotes the true value of 9, the theory behind the 

GMM states that 9^ ^ is approximately distributed as Normal(0, Vt/T) with 



Vn 



d_ 1 
M \ T 



Sir, 



8=0 K 



-i d I 1 
dO 



6=0 



(i) 



(fexr) 
it 



(rxk) 



For the case under consideration we have 9 t = (/i, Ty , y^ ,yo,ay,Tz, ,zq,<j%, pxy ,Pxz, Pyz) , 
while the orthogonality conditions can be obtained computing the lowest order mo- 
ments of returns, the leverage correlation, and the squared return autocorrelation 

AX t 



E[h(0,W t )]=E 



I AX, 



AX t AX f 2 +Ai - At 2 (C 2 ,o + 2Ci A + C , 2 ) 2 £{L'At; t) 

AX t AX? +LAt - At 2 (C 2 , + 2d.! + C , 2 ) 2 C(L"At; t) 
AX?AX 2 +K , At - At 2 x r.h.s. of Equation (C.7) for r = K' At 

AX 2 AXf +K „ At - At 2 x r.h.s. of Equation (C.7) for r = K" At 



where AX t — \x\St+At — hi St — pAt, At — 1/250 yr. In the previous equation, the 
conditions referring to the return- volatility and the squared return correlations depend 
on the four positive integers L' < L" , K' < K" . The choice of these values can be 
made based on a prior analysis of the time scales of the correlation functions, and 
will be detailed later on. Thus, the dimension r of the vector h.(9,W t ) reduces to 
4 + L" — L' + 1 + K" — K' + 1. From an econometric point of view the problem of 
the estimation of parameters is cast into a sound statistical framework. By means of 
GMM we can obtain an estimate of central values and associated statistical uncertainty 
for all the unknowns of the problem. However, the quantity to be optimized is 
highly non linear, the optimization procedure of the twelve dimensional problem 
is per se problematic, and finding a solution under blind search can be extremely 
demanding. For this reason, we prefer to proceed by invoking some reasonable 
arguments concerning the nature of the problem under study. The starting point of our 



§ When the process h(0, Wt) for t = 1, . . . , T is serially correlated, the Newey-West estimate for SIt 
can be used, please refer to equation 14.1.19 in 33 for further details. 
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heuristic is the observation that, until now, we have devoted little attention to the role 
played by the parameter t . In principle it could be treated as an unknown parameter, 
however its role is quite different from that played by the others. Since it mainly 
determines the regime of the factors processes, we assume to —oo as done in the 
previous work 16J. Said differently, we assume that the data we are observing reflect 
stationary realizations of Y t and Z t . Under this regime, mean-reverting processes do 
not depend on the initial time values j/o and zq any more, and we identify j/oo with 
2/0) and z^ with z . Moreover, both Y t and Z t are unobserved processes reflecting 
the presence in the market of investment strategies with heterogeneous time horizons. 
Even though this assumption could be relaxed, it is plausible to assume that the 
Brownian motions driving those processes are uncorrelated. If we fix pyz — the 
problem greatly simplifies since all F mi „ reduce to + F„, and all terms C mj „ 
split into C m fl x Co >n . In [TH] we prove that pxy < suffices in order to reproduce 
the leverage effect. Since we do not want that the incorporation of the extra factor 
Z t has a relevant impact on the leverage, we fix pxz equal to zero. Finally, the 
considerations that follow Equation ([5| in Section [2] have clarified the way the tail 
exponent of the distribution of the volatility factors is responsible for the divergence 
of the moments of X t . If Y t and Z t were characterized by two different tail exponents, 
the order of the first divergent moment of X t should be determined by the lowest of 
them. In this respect the role played by the highest exponent would be spoiled by 
the other one. We therefore assume that the stationary distributions of Y t and Z t 
have the same shape parameter v = vy =Vz- Now the reduced vector of parameters 
reads t = (jj,, y^, Zoo, Ty, tz, Pxy, v), while the orthogonality relations simplify. For 
instance, the first four relations reduce to 



E [AX t ] - 
AX t \ 



0, 



E 



E 



E 



2At 



= 0, 



(AX t ) -(y 00 + z co yAt + 



v — 2 



\AX\ 



8At 3 (v - If 
7T 0-3)0-2) 



E 



^^SAt^v — 1 

7T V — 2 



iVooZo 



= 0. 



and the numerator of the leverage for positive r becomes 

1 



PXY] 



T Y {y - 1) 



2,0 L/ 0,1 qy°° O l,0 o 0,l 



exp 



- 1 



(9) 



v- 2 



v 



1 



[C^oCo^l + ^1,0^0*2 — z oo (C^O + Cl'oCo!l)] ex P 

1 



TZ) T C 



{jZ-\y°o + ^oc) (C^o + Ci!oC^) exp ( 



where the superscript st stands for the stationary regime corresponding to to — > — oo, 
and we recall that T£ = (y — l)ry j{y — 2). Even though the leverage correlation 
introduces a superposition of three exponential functions, we have seen at the end 



of the Section 3.1 that the characteristic exponents are of the same magnitude and 
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are all dominated by T£. For this reason, and recalling that the typical decay time 
for the leverage is smaller than one hundred days, to perform the optimization we fix 
L' = 1, L" = 250, and solve the problem for the first 254 orthogonal relations. This 
greatly enhances the convergence of the numerical algorithms. Once an estimate of 
T£ is found, we fix K' equal to two times the integer part of re, K' = [2tc\ [jj] and 
K" = 250, and we perform the final optimization on the entire set of 254 + K" — K' + 1 
orthogonal relations. The latter relations, which correspond to the volatility structure, 
should be fixed using the whole range for the lagged correlation, and not just the lags 
above the leverage. However, we preliminarily perform the numerical optimization 
with K' running from one to the mentioned level and we see that the fit worsens. The 
reason comes from the very structure of the expression which appears in the r.h.s. of 
Equation (C.7). As commented in Section 3.2 the lagged correlation is dominated 
by two time scales, the first one is of order of T£, while the second one is the long 
run component. When K' = 1, we are fitting the behaviour of the whole curve, but 
for low lags it is largely determined by the scaling of the leverage. The results of 
our numerical explorations show that the low lags part of the curve is hard to be 
reproduced while the long run component is always well described. When we reduce 
K' to one the optimizer tries to catch up at short lags but at the cost of an even worst 
distortion of the autocorrelation for intermediate lag values, please refer to Figure [2] 
for a visual comparison of the cases K' = 1 and K' = [2tc\ . In light of these results 
we decide to perform the GMM favouring the long run behaviour. An alternative 
and less abrupt approach could be a weighted optimization on the entire curve with 
weights in the matrix fir- that change smoothly from zero to one from low to high lag 

values. Finally we perform the GMM just one time and we obtain 6 ^ and Vt/T. 
In order to compute a consistent estimate of o~ Y and the associated confidence level, 

we can extract a random sample from Normal(^ , Vt/T) and obtain a statistics of 
a Y through the relation 2/(ty{v — 1)). We proceed in an analogous way for o\, and 
for T£ = Tyiy — 1)/ {v — 2). The time series on which we perform the analysis is the 
same used in [16] , and it consists of a data set from the Standard & Poors 500 index 
daily returns from 1970 to 2010. This allows to evaluate the ability of the extended 
model to capture the persistence of the volatility, not only in absolute terms but also 
in comparison with the previous estimate from a simpler model. In Table [T] we report 

the central values the standard errors &t = \J diag(VT/? 1 ), and the correlation 

structure p T = ~Vt/(To , t^'t) f° r au the parameters. As far as the other relevant 
parameters of the model are concerned, we have a Y = 9.9± 1.9, a\ = 1.58±0.07, and 
T£ = 0.10 ± 0.02 yr. The new values confirm the goodness of the estimate provided 
in [16], in particular the value of pxy is strictly negative and the level of the tail 
exponent v predicts the divergence of moments higher than the fourth one. More 
interesting to comment is the relationship between the different time scales involved 
in our process. Indeed, the shortest time scale corresponds to the typical relaxation 
time of Y t , which is found to be equal to 0.07 ± 0.01 yr and is therefore dominated 
by the leverage time scale 0.10 ± 0.02 yr (to be compared with the old estimate for 
T£ in [16 which is 0.09 yr). The new time scale Tz for the process Z t is found to be a 
factor of six larger than that of Y t . In Figures [l] and [2] we plot the leverage function 

|| In the following we adopt standard mathematical notation [-J for the integer part function. 
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Figure 2. Empirical volatility autocorrelation function of the daily returns of 
the S&P500 index 1970-2010 (data points), and analytical descriptions: bold and 
dotted lines, new expressions with GMM estimates for different values of K'\ 
dashed line, formula and values of parameters as in I16| . 




T 



and the normalized autocorrelation of squared returns. The exponential decay of the 
leverage is described correctly by the analytical formula, and no relevant differences are 
noticeable with respect to the description obtained via the model introduced in |16) . 
Different considerations apply to the persistence of the volatility as predicted by the 
extended model. The presence of the slow volatility factor Z t introduces a longer 
time scale allowing to capture the long range memory of the autocorrelation function. 
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Table 1. Estimated values of the parameters (K' = [2t£J) from daily returns of 
the S&P500 index 1970-2010 . 
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This is evident from the comparison between the dashed line, corresponding to the old 
model, and the bold one, corresponding to model Q. Our results demonstrate the 
ability of a multi-factor approach to stochastic volatility to effectively describe several 
phenomena. In particular, even in the simplest version of a two factor model, it is 
able to capture the emergence of multiple time scales for the volatility autocorrelation 
as well as the exponential decay of the return-volatility correlation. The measured 
value for the tail parameter v is coherent with the internal consistency of the model 
requiring v to be greater than four (in order for Equation (C.7) to converge in the 
stationary limit). In particular, v — 4.15 predicts an hyperbolic decay of the daily 
return distribution which captures correctly the non Gaussian probability of extreme 
events in the real data. 



5. Conclusions 

In this work the model for the description of financial stylized facts proposed in [16] 
is amended from the unrealistic fast decay of the volatility autocorrelation. This is 
achieved introducing an extra stochastic factor driving the volatility. In principle the 
number of factors could be increased at will, but the analytical tractability of the 
resulting model would be hardly exploitable. The intuition behind this generalization 
traces back to the early empirical analysis of the FX market in [9] and the model in [10] , 
where the role played by heterogeneous investors is strongly emphasized. Evidences 
from these papers are rooted in the econometric analysis of publicly available financial 
time series, but a convincing micro-founded model is still lacking. Access to electronic 
order book data and to agents' identifiers would allow to estimate the individual 
components of this heterogeneity. An even approximate estimation of the distribution 
of typical investment horizons from this information would provide a valuable trader- 
based foundation. 

With respect to previous approaches and analyses of continuous-time stochastic 
volatility models, we believe that the calibration procedure proposed here represents a 
further improvement and fulfils the desirable requirements of statistical soundness. At 
the same time, it also allows to focus on those facts which are established as relevant 
for the description of financial data. In particular, we pursue an heuristic approach to 
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optimization that, reducing the dimensionality of the parameters space, retains only 
those ingredients which are actually needed to capture the aforementioned empirical 
evidences. 

The stochastic volatility models discussed in this paper and in |16j are inspired 
exclusively by the quest for a realistic description of financial data. In this quest, we 
focus on continuous time modelling and do not consider aspects relating to possible 
applications in the financial sector. In this respect, we should mention that important 
progresses have been obtained by discrete time models. In particular, models with 
multiscale ARCH volatility can also accommodate many stylized facts, like the fat 
tails of returns, and reproduce consistently the dynamics of realized volatility, also 
delivering accurate forecasts of the latter [34]. They also prove to be very flexible 
tools for efficient option pricing and hedging purposes [35] . 

In our case, the emergence of power-law tails in the return distribution complies 
with past empirical analysis |36j . but also poses serious limitations to the usage 
of the model in the context of option pricing. This is certainly true for vanilla 
instruments, whose payoff grows exponentially with the log-price. On the other hand, 
our framework could deliver, in principle, better estimates of the role played by rare 
events for market risk evaluation. 

On the whole, we believe our model achieves a remarkable degree of realism, higher 
than previous attempts in continuous-time stochastic volatility modelling, yet allowing 
for important analytical derivations, e.g. that of the moments of the return probability 
distribution. Fairly enough, this comes at the price of elaborate manipulations and 
Monte Carlo simulation would still be due for most financially relevant applications. 
Nonetheless recent advances in fast computing, e.g. GPU based numerical techniques, 
could offer a promising scenario in this regard. 
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Appendix A. Probability density function of price returns: moments 
computation 

Application of Ito's Lemma to the function X\ readily provides 



^ Jo J 

and the same Lemma proves that the correlation functions between integer powers of 
Xt, Yt, and Z t satisfy the following differential equation 




(A.l) 



^E [X l t Y™Z?] = F m ,„ E [X l t Y t m Z?] + A Y m E [X^" 1 ^] 

+A z n E [X^zr 1 ] + - 1)E [x\- 2 Y t m Z? (Y t + Z t f 
+ l(mp XY a Y + np xz a z ) E [X^Y^Z? (Y t + Z t )] , 



(A.2) 
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where the constants F m ^ n are defined right after Equation ([5j) and Aj n = mKyyoc 
and = nn z z ao . Previous equations correspond to a system of nested linear 
ordinary differential equation (ODE), which can be solved recursively starting from 
the lowest order of I, m, and n, and whose solution involves integration of the two 
point correlations C m . n (t;to) = E[Y t m Z™] ^ From application of Ito's Lemma we 
obtain 

d (Y rn Z n ) = [F%Y m Z n + F^Y m Z n ] dt 
+ [A Y l Y m - 1 Z n + AlY m Z n - l \ dt + p YZ mn^a Y a 2 z Y m Z n dt 

Y i „ . /„2 \rm rznjTjrZ . 



-m\/a Y Y m Z n dW t 



n\/a 2 z Y m Z n dW t k , 



E [Y™Z?] = F* + F* + p YZ mnJa Y a% E [Y t m Z?] 



taking expectation, and differentiating w.r.t time we derive the following ODE 
d. 
dt" 

+AlE [Y^-'Z?] +A*E [Y^zr 1 ] 
For instance, for the case m = n = 1 we have 
d. 



E [Y t Z t ] = (F* + Ft + p YZ Ja Y a% E [Y t Z t ] 



dt 

+A Y E [Z t ] + Af E [Y t ] , 
where the mean values read 

E[Y t ] 

with to < the starting time of the factors processes. More generally, by iterative 
solution it can be verified that C mj „ admits the expansion 



A\ 
FY' 


|_ e ^i i '(t-to) 


yo + 


AY 
FY 


Af 

FY' 


i_ e *f(*-*o) 


zo + 


Af 
FY 



C m , n = E [Y™Z?\ = J2 E kf^e F ^-^ , (A.3) 

z=0 j=0 

where the coefficients depend on the initial conditions yo, z$ and satisfy the recursive 
relations which follow 



, (m,n) 
i<m : j <n 



k 



(m,n) 



<m,n 



m,n) 
m,n 



AY u(.m-l,n) , , Z Am,n-\) 
m i,j + n K i,j 

F, — F- ■ 

± m,n ± i,j 
AY T.{m-l,n) 



F„ 



F,,_ 



n m j 
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± {, 



F„ 



m — 1 n 



,(m-l,n) 



m ri— l ,(m,n— 1) 



^EE 



i=0 j=0 



Fi. 



1 In the following we drop the dependence on t and tg. 



F, 



(A.4) 
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We notice that t he m oments Hm(t) = E [Y t m ] and n%(t) = E [Z™] are specific cases 
of the expansion (A.3), whose coefficients are given by the column vector (&j™'°^)i< m 
and the row vector (fco°/^).7<n; while in general the set of coefficients k^' n ^ can be 
cast in a (m + 1) x (n + 1) real matrix. For instance, for the case 62,1 = E [Y 2 Z t ~\ we 
obtain 

1 



1.(2.1) 
^0,0 



1.(2,1) 
K 0,l 



F2.I 



A Y k {1 > 1) + A z k {2 ' 0) 

A 2 ^0,0 + A l ft ( 



^0,0 



k 



(2,1) 
1,0 

1.(2.1) 
"-1,1 

k^ 

h 



F2.1 — F x z 
1 

^2,1 — f y 

A Y k ihl) 

1\ 2 K-y -y 



F2.1 — F] 



A Y k {1 < 1) 



Sly «,]_ o 



1,1 



2.0 

(24) 



A z k {2fi) 
F2.1 — F 2 



v 2,l 



1.(2,1) 
v 0,0 



-,(2,1) 
c 0,l 



u(2,l) 
"1,0 



u(2,l) 
v l,l 



+ k 



(2,1) 
2,0 



Ultimately, given the expansion (A.3) and by inspection of ( |A.l 1, we recognize that 
the moments of X can be cast in the form of Equation ([5| . 



Appendix B. Computation of the return- volatility correlation 

The results which follow are derived under the assumption that pxz is equal to zero. 
The numerator of the return-volatility correlation £(r; t) = E [dX t dA t 2 +T ] /E [dA t 2 ] 
can be cast in the form 

E [dX t dXf +T ] = E [(Y t + Z t ) (Y t+T + Z t+r f Cf ] dt 2 . 

Here, adopting the same convention of [25], we formally [^] express the Wiener 
increment as dW^ = Ct dt, where C,^ is a Gaussian noise with zero mean and variance 
E [(C* X ) 2 ] = 1/dt. Novikov's theorem [37] [38 allows to compute the expectation 
involving £^ , giving us 

E [dX t dX 2 +T ] 



dt 2 
E 



%Pxy \/^y^( t ) exp (— Kyr) x 



[Y t (Y t + Z t )(Y t+T + Z t+T )] exp 



o$A t W x (r) 



(B.l) 



where we define A t M / (r) = J f dW s . We refer the interested reader to Section IV 
in |16| for further details regarding the derivation of the previous equation. The right 
hand side of (B.l ) can be split into four pieces proportional to the expectations 



f Y YY(T,t)=E 



Y 2 Y t+T exp 



a 2 Y A t W Y (r) 



+ This convention is somewhat unusual, but admissible as the distributional assumption on ti* 
guarantees that dWt is still Gaussian with zero mean and variance equal to df . This choice proves to 
be convenient to develop the following calculations. 
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/yy Z (r,i)=E 

fYZY(T,t) =E 

fyzz(r,t) =E 



Y t Z t+T exp 
Yt Z t Y t+T exp 
Yt Z t Z t+T exp 



a y A t T^ Y (r) 
^ A t W y (r) 
'4 A t W Y (r) 



Following the approach discussed in Appendix B of [16] . it is possible to show that 
they satisfy the relations 



fYYY(r,t) - (cry - k y ) Jyyy (t', t) exp 



-(r-r') 



dr' 



exp ( -y-r ) [C 3j o + KyyooTC 2 ,o] 



Jyyz{t, t) + Kz / ,/Vyz(t', t) exp 



-(r-r') 



dr' 



exp ( ~|~ r ) [C2,oC ,i + kzZootC 2i0 ] , 
fYZY{r,t) - (cr y - Ky) / /y Z y(T"',i)exp 



-(r-r') 



dr' 



exp ( -y-r ) [C 2 , Co,i + KyyooTC li0 Co,i] 



fYzz(T,t) + K Z Jyzz{t' ,i)exp 
Jo 



(r-r') 



dr' 



CXp ( ^yT J [Ci )0 C ,2 + KzZaoTC-LflC^i] , 

corresponding to a set of Volterra integro-differential equations of the second kind. 
Their solutions are known in closed- form, and after plugging them in Equation (B.l), 
the final expression of the leverage correlation reads 



2 Px y^H{t) 



(C2,o + 2Ci,oCo,i + Co, 2) 

KYUoo 



C3,o + C2 .qCo 1 + 



(C2,0 + Ci,oCo,l) 



a Y — Ky 

+ [C 2 ,oCo,i + Ci, 2 - Zoo (C 2 , + C 1)0 C 0il )] exp 

KYVoo 



exp 
2 



y 



3c? 



Ky - K Z T 



Ky T 



0~ Y ~ Ky 



Zoo I (C2,0 + Ci : qCo,i) 



exp 



Ky T 



(B.2) 



A meaningful comparison of the previous expression with real data requires to take 
its stationary limit for t — > —00, which amounts to replacing C m ,oCo,„ with the 
asymptotic values C^ Cg* n . 
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Appendix C. Computation of the squared return correlation 



Resorting to the same parametrization of the Wiener variation adopted in [Appendix] 
|B]we have 

E [dX? dXf +T ] = dt 2 E \(Y t + Z t f (Y t+T + Z t+T f dW t x Cf 1 



dt 2 E 



(Y t + Z t f (Yt+r + Z t+T ) 2 +0(dt 3 ). (C.l) 



In order to compute the autocorrelation function of squared returns, the quantities 
/ t (m '"' M) (r) = E [Y t m Z?Y? +T Z% +T ] indicating the r-lagged correlation have to be 
evaluated. The relevant cases correspond to p, q < 2, and below we detail the 
corresponding exact results, all of which are obtained replacing the process Yf, T Z^, T 
with its integral representation from time t to time t + r. 



Computation of f\ 



(m, n,l,0) 



(t) = E [Y t m Z?Y t+T } . It is readily verified that /, 



(m, n,l,0) 



(t) 



is solution of a linear ODE, giving 



ft 



(m,n, 1,0) 



(r) 



(C.2) 



Computation of f( m ' n >°>V ( T ) = g [y™Z t "Z f+T ] . In much the same way we have 

A? 



fi 



(m,n, 0,1) 



A z 

( T ) = T^Z^ m > n 
^1 



Computation of f^ 2 ^ ( T ) = E [r™Z t "y t 2 +r ] . After replacement of Y 2 +T1 we can 



write 



/■(m,n, 2,0) / \ ^(m,n,2,0) 



(0) 



(m,n, 2,0) 



(T')dT' 







p(m,n, 1,0) / / 



(r')dr'; 



further, we can replace the solution (C.2 1 for /, 



,1,0) 



(V) in the second integral, 



leading to straightforward integrations of exponential functions of r. Finally, we are 
left with 

A\A\ 



fi 



,2,0) 



(r) 



F 2 F 1 



Cm+2,n + 



^m+l,n "T ^m^n 



r 2 r l 



4> 

r i_ 1 r 

b 2 



Y 



(C.3) 



(m,ra,0,2) 



(t) = E [Y t m Z\ l Z 2 +T ] . As before, after replacement of the 



Computation of / ( 

parameters for the dynamics of the Z 2 , T process, we get to 



fi 



(m,n, 0,2) 



F 2 Z F? 
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2 * 1 

C m , rl +2 + 



Af, 



pZ _ pZ I Cm,n+l + p Z C m>n 



Fi - F? 



— ( r a- A ^ r 



(C.4) 



Computation of ^ m ' n,1 ' 1 ^ 7 -) 
process F t Z t is given by 

d(F t Z t ) = (F 14 F t Z t + Af Y t + A\Z t ) dt 



E[Y t m Zl l Y t+T Z t+ r}- The evolution of the joint 



u\Y t Z t AWj 
and substitution inside the expectation gives 



n 



(m,n, 1,1) 



(r) 
Af 



\A\Af 













*M - ^ i 



if 



m+l,n 



fy 



if i - if 



r 1 x ^Lr 



Jl.1T 



AfAf 



a 



Af 



m+l,n+l T y 

2if i - if - if 



Cm+l,n + 



AY 



a 



m,n+l 



(C.5) 



^1,1(^,1 -*iO(*U-*i z ), 

As expected from the structure of model fil), and as confirmed by all previous 



examples, it is clear that the functions f( m ' n ' p ' q> ( T ) admit a general expansion reading 



/ (-^)( T)= ^^ /l (-.™ (t)e ^ 
i=i j=i 

(i) can be computed exactly. Coming back to Equation 



^ 



where the terms fi^ m ' n ' p < q > 
(C.l ) we have 



E[dX?dX? +T ] 
dt 2 



/t (2,0, 2l 0) (T) + /t (0,2,2,0) (r) + 2/ (0,2,2,0) (r) + ^(2,0,0,2)^ + ^(0,2,0,2)^ 

+2 / t (0 < 2 <°< 2) (r) + 2 [ jf (r) + / t (0 < 2 < M) (r) + 2 / t (0 ' 2 ' M) (r)l . (C.6) 



By means of Equations (C.3)-(C.5), and after defining the auxiliary variables 
T\ = G'2,0 + Co, 2 + 2(71,1 1 ^2 = £'3,0 + C1.2 + 2C2,i , 

T 2 = Co, 3 + C2,l + 2Ci,2 , T 3 = C4,o + C2,2 + 2C34 ' 
^3 = Co, 4 + C2,2 + 2Ci,3 i T4 = C3,i + Ci,3 + 2C2,2 ; 

we can write the following final expression 
E [dX? dX 



dt 2 



pr pr ■ pZpZ 
e 2 r l r 2 r l 



Fi, 







(w- 





Ti 



A Y 



Af 



F*-F? ifi-if 
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To* 



T 3 



Ff 



F z 



Ff 



if 



2e 



-Fl.lT 



T 4 

2F M 



F/ -F/ 
Ff-Ff 



T 2 + ^Tx 



Fl! - Ff 



T* 



rT* 



To 



FY 



Ff 



F 1 , 1 (F 1)1 -Ff)(F 1)1 -Ff)J • (CJ) 

Ultimately, evaluation of the volatility autocorrelation ([6]) requires to compute 
Var[dX 2 ] = E [dX 4 ] - E [dX 2 ] 2 which is given by 

3 (C 4 . + 4C 3 ,i + 6C 2 , 2 + 4Ci, 3 + C a) dt 2 - (C 2)0 + 2C M + C , 2 ) 2 di 2 
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